If a periodic system contains charges interacting via the \(\frac{1}{r}\) Coulomb potential, direct summation of the interactions
(1)\[E = \sum_{(i,j)} \frac{1}{4\pi\epsilon_0}\frac{q_i q_j}{r_{ij}},\]
where the sum is over pairs of charges \(q_i, q_j\) (charges of the entire system, not just the simulation cell) and the distance between the charges is \(r_{ij} = |\mathbf{r}_j - \mathbf{r}_i|\), does not work in general because the sum (1) converges very slowly (it actually converges only conditionally). Therefore truncating the sum may lead to severe errors. More advanced techniques must be used in order to accurately evaluate such sums.
This class represents the algorithms used for evaluating the \(1/r\) sums. It wraps the summation parameters and activates the summation of Coulomb interactions. If an instance of CoulombSummation is given to the Pysic calculator, Coulomb interactions between all charged atoms are automatically included in the calculations, regardless of possible Potential potentials the calculator may also contain. Otherwise the charges do not directly interact. This is due to two reasons: First, the direct Coulomb interaction is usually always required and it is convenient that it is easily enabled. Second, the specific potentials described by Potential are evaluated by direct summation and so the Coulomb summation is separate also on algorithm level in the core.
Sometimes, you may want to scale the effective charges before calculating the Coulomb sum. Especially, you may want to exclude some atoms from the long range summation. This can be done by giving the CoulombSummation a list of scaling values, one per atom. The actual charges of the atoms are then multiplied by the given scaling values before the Coulomb potential is calculated. If a scaling value is 0, the corresponding atom is always treated as if it had no charge. Note though, that scaling with unequal scaling constants may lead to the cell being effectively charged.
Using the Python map function is a convenient way to generate such atom-by-atom lists. For instance, if you want to generate a list by element:
>>> atoms = ase.Atoms('H2O')
>>> def charge_by_elem(elem):
... if elem == 'H':
... return 0.1
... elif elem == 'O':
... return -0.2
... else:
... return 0.0
...
>>> system.set_initial_charges(map(charge_by_elem, system.get_chemical_symbols()))
>>> charges
[0.1, 0.1, -0.2]
As the summation algorithms are parameter dependent, one should always check numeric convergence before real simulations. As a first guess, the utility function pysic.interactions.coulomb.estimate_ewald_parameters() can be used for estimating the parameters of the Ewald method.
Below is a list of summation algorithms currently implemented.
The standard technique for overcoming the problem of summing long ranged periodic potentials is the so called Ewald summation method. The idea is to split the long ranged and singular Coulomb potential to a short ranged singular and long ranged smooth parts, and calculate the long ranged part in reciprocal space via Fourier transformations. This is possible (for a smooth potential) since the system is periodic and the same supercell repeats infinitely in all directions. In practice the calculation can be done by adding (and subtracting) Gaussian charge densities over the point charges to screen the potential in real space. That is, the original charge density \(\rho(\mathbf{r}) = \sum_i q_i \delta(\mathbf{r} - \mathbf{r}_i)\) is split by
\[\begin{eqnarray} \rho(\mathbf{r}) & = & \rho_s(\mathbf{r}) + \rho_l(\mathbf{r}) \\ \rho_s(\mathbf{r}) & = & \sum_i \left[ q_i \delta(\mathbf{r} - \mathbf{r}_i) - q_i G_\sigma(\mathbf{r} - \mathbf{r}_i) \right] \\ \rho_l(\mathbf{r}) & = & \sum_i q_i G_\sigma(\mathbf{r} - \mathbf{r}_i) \\ G_\sigma(\mathbf{r}) & = & \frac{1}{(2 \pi \sigma^2)^{3/2}} \exp\left( -\frac{|\mathbf{r}|^2}{2 \sigma^2} \right) \end{eqnarray}\]
Here \(\rho_l\) generates a long range interaction since at large distances the Gaussian densities \(G_\sigma\) appear the same as point charges (\(\lim_{\sigma/r \to 0} G_\sigma(\mathbf{r}) = \delta(\mathbf{r})\)). Since the charge density is smooth, so will be the potential it creates. The density \(\rho_s\) exhibits short ranged interactions for the same reason: At distances longer than the width of the Gaussians the point charges are screened by the Gaussians which exactly cancel them (\(\lim_{\sigma/r \to 0} \delta(\mathbf{r}) - G_\sigma(\mathbf{r}) = 0\)).
The short ranged interactions are directly calculated in real space
\[\begin{eqnarray} E_s & = & \frac{1}{4 \pi \varepsilon_0} \int \frac{\rho_s(\mathbf{r}) \rho_s(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} \mathrm{d}^3 r \mathrm{d}^3 r' \\ & = & \frac{1}{4 \pi \varepsilon_0} \sum_{(i,j)} \frac{q_i q_j}{r_{ij}} \mathrm{erfc} \left( \frac{r_{ij}}{\sigma \sqrt{2}} \right) - \frac{1}{4 \pi \varepsilon_0} \frac{1}{\sqrt{2 \pi} \sigma} \sum_i^N q_i^2. \end{eqnarray}\]
The complementary error function \(\mathrm{erfc}(r) = 1 - \mathrm{erf}(r) = 1 - \frac{2}{\sqrt{\pi}} \int_0^r e^{-t^2/2} \mathrm{d}t\) makes the sum converge rapidly as \(\frac{r_{ij}}{\sigma} \to \infty\). The latter sum is the self energy of each point charge in the potential of the particular Gaussian that screens the charge, and the sum runs over all charges in the supercell spanning the periodic system. (The self energy is cancelled by the long range part of the energy.)
The long ranged interaction
can be calculated in reciprocal space by Fourier transformation. The result is
\[\begin{eqnarray} E_l & = & \frac{1}{2 V \varepsilon_0} \sum_{\mathbf{k} \ne 0} \frac{e^{-\sigma^2 k^2 / 2}}{k^2} |S(\mathbf{k})|^2 \\ S(\mathbf{k}) & = & \sum_i^N q_i e^{\mathrm{i} \mathbf{k} \cdot \mathbf{r}_i} \end{eqnarray}\]
The sum in \(E_l\) runs over the reciprocal lattice \(\mathbf{k} = k_1 \mathbf{b}_1 + k_2 \mathbf{b}_2 + k_3 \mathbf{b}_3\) where \(\mathbf{b}_i\) are the vectors spanning the reciprocal cell (\([\mathbf{b}_1 \mathbf{b}_2 \mathbf{b}_3] = ([\mathbf{v}_1 \mathbf{v}_2 \mathbf{v}_3]^{-1})^T\) where \(\mathbf{v}_i\) are the real space cell vectors). Likewise the sum in the structure factor \(S(\mathbf{k})\) runs over all charges in the supercell.
The total energy is then the sum of the short and long range energies
\[E = E_s + E_l.\]
If the system carries a net charge, the total Coulomb potential of the infinite periodic system is infinite. Excess charge can be neutralized by imposing a uniform background charge of opposite sign, which results in the correction term
This correction is applied automatically.
Forces are obtained as the gradient of the total energy. For atom \(\alpha\), the force is
(There is no contribution from \(E_c\).) The short ranged interactions are easily calculated in real space
where \(\hat{r}_{\alpha j} = \mathbf{r}_{\alpha j} / r_{\alpha j}\) is the unit vector pointing from atom \(\alpha\) to \(j\). The long range forces are obtained by differentiating the structure factor
Below is a list of methods in CoulombSummation, grouped according to the type of functionality.
Class for representing a collection of parameters for evaluating Coulomb potentials.
Summing \(1/r\) potentials in periodic systems requires more advanced techniques than just direct summation of pair interactions. The starndard method for evaluating these kinds of potentials is through Ewald summation, where the long range part of the potential is evaluated in reciprocal space.
Instances of this class are used for wrapping the parameters controlling the summations. Passing such an instance to the Pysic calculator activates the evaluation of Coulomb interactions.
Currently, only Ewald summation is available as a calculation method.
Parameters:
Creates a dictionary of parameters and initializes all values to 0.0.
Sets a given parameter to the desired value.
Parameters:
Sets the numeric values for all parameters.
Parameters:
Sets the numeric values for all parameters.
Equivalent to set_parameter_values()
Parameters:
Set the list of scaling parameters for atomic charges.
Parameters:
Sets the summation method.
The method also creates a dictionary of parameters initialized to 0.0 by invoking initialize_parameters().
Parameters:
Names of the summation methods. These are keywords used for setting up the summation algorithms.
Short descriptions of the parameters of the summation algorithm.
Names of the parameters of the summation algorithm.
Returns a tuple containing a good initial guess for Ewald parameters in the order real_cutoff, k_cutoff, sigma, epsilon.
The returned values are (real_cutoff, k_cutoff, sigma, epsilon). Epsilon is always 0.00552635 \(\varepsilon_0\) in units of \(\frac{e^2}{eV A}\) Real cutoff can be given by the user and should be short enough to make the real space summation efficient (note that this affects neighborlisting and thus also other real space summations). Sigma and k_cutoff are determined by simple scaling rules where \(\sigma \sim r_\mathrm{cut}\) and \(k_\mathrm{cut} \sim \sigma^{-1}\). The scaling constants are determined by the accuracy setting.
Note that the given parameters are not analysed in any way - they are only a first guess. You should always test the parameters for accuracy and speed before production simulations.
Parameters: